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Abstract 

In this paper we present a computation of the mean first-passage times both for a random walk 
in a discrete bounded lattice, between a starting site and a target site, and for a Brownian motion 
in a bounded domain, where the target is a sphere. In both cases, we also discuss the case of two 
targets, including splitting probabilities, and conditional mean first-passage times. In addition, we 
study the higher-order moments and the full distribution of the first-passage time. These results 
significantly extend our earlier contribution [Phys. Rev. Lett. 95, 260601]. 

PACS numbers: 05.40.Jc,05.40.Fb 
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I. INTRODUCTION 



The time it takes to a random walker to go from a starting site to a target site, the so 
called first-passage time (FPT), is an especially important quantity that underlies a wide 
range of physical pro cesses Q, [4]. Indeed, numerous real situations, such as diffusion limited 
reactions Q or animals searching for food {4], can be rephrased as first-passage problems. In 
all these situations, the FPT is a limiting factor. As a consequence, it is crucial to determine 
how this quantity depends on the parameters of the problem. 

Among these parameters, geometrical factors turn out to be determining. For example, 
the mean first-passage time (MFPT) between a starting site and a target site for a 2D 
random walker is infinite if the walk is not bounded. On the contrary, it becomes finite as 
soon as the walk is confined. But how does the MFPT depend on the confining surface? 
In fact, the answer to this general question appears as a difficult task, because explicit 
determinations of FPT are most of the time limited to very artificial geometries, such as ID 
and spherically symmetric problems j^j]. 

However, in most of the real situations, the searcher performs a random walk in more 
general confining geometries. This is for example the case in biology, where biomolecules 
often follow a complicated series of transformations, which are located at precise parts of 
the cell. Determining the influence of the shape of the cell on the FPT actually appears as 
a first step in the understanding of the global kinetics of the process. 

This question of determining first passage properties in gene ral confined geometries has 
raised a growing attention (see for example 

important results have notably been obtained. First, in the case of discrete random walks, 
an expression for the mean first passage time (MFPT) between two nodes of a general 
network has been found [16( . However, no quantitative estimation of the MFPT was derived 
in this paper. Second, the leading behavior of MFPT of a continuous Brownian motion at 
a small absorbing window of a general reflecting bounded domain has been given 17|, Il8 |. 
These studies have even been extended to a situation with a deep potential well, leading 
to a generalization of the Kramers formula 19) . In the case when this window is a small 
sphere within the domain, the behavior of MFPT has also been derived 2(|. This result is 
rigorous, but does not give access to the dependence of the MFPT with the starting site. 
Very recently 2JJ, we have proposed a novel approach which allowed us to propose 
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accurate estimations of first passage times of discrete random walks in confined geometry. 
Preliminary results concerning a continuous Brownian motion have also been announced. 
The main purpose of this paper is to provide a detailed analysis of this continuous case, 
relevant to many real physical situations. In addition, we extend our previous work in 
several directions, for both discrete and continuous cases: the complete distribution of FPTs 
is obtained; extra quantities, as conditional MFPT in the case of several targets or mean 
exit times by a small aperture of a general reflecting bounded domain, are derived. 

The paper is structured as follows. In Sec. [Til we first present the computation method of 
FPTs in the case of random walks on discrete lattices. This study includes the obtention of 
the MFPT, a comprehensive derivation of the expression of the higher order moments as well 
as the complete distribution of the FPT, whose physical meaning is extensively analysed. 
The situation with two competitive targets is also studied, and we compute MFPT, splitting 
probabilities, and conditional MFPT. 

In Sec. IHIl we extend all these results to the case of a continuous Brownian motion, and 
detail the specific difficulties encountered in this case. 

The explicit results obtained in Sees. HIJ and IHIl involve pseudo-Green functions of a 
Laplace type operator, with given boundary conditions. The Appendix [X] is devoted to the 
evaluation of these pseudo-Green functions. For several domain shapes, an exact formula 
can be obtained, which gives, for the quantities computed in the article, exact explicit 
expressions in the discrete case, or accurate approximations in the continuous case. For 
other domain shapes, basic approximations are proposed. 

These results are briefly summarized in Sec. IIVI with a discussion of the important 
parameters to take into account and of the qualitative behavior of the MFPT 



II. RANDOM WALKS ON DISCRETE LATTICES 



A. Mean first-passage time 

Let us consider a point performing a random walk on an arbitrary bounded lattice with 
reflecting boundaries. We want to compute the MFPT (T) of the random walker at target 
site T, starting from a site S at time 0. We summarized this computation in a previous 



article 21]. 
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However, since it is the basis of all the developments explained in this article, we found 
useful to give it here in full detail, with the addition of several necessary precisions. 

Our method is based on a formula given by Kac 22|, concerning irreducible graphs, such 
that any point can be reached from any other point. An irreducible graph admits a unique 
stationary probability 7r(r) to be at site r (physically, this is the probability for a particle 
which has been in the domain for a long time to be at site r. If the transition probabilities 
are symmetric this stationary probability is uniform.). We consider random walks starting 
from an arbitrary point of a subset £ of the lattice, chosen with probability 7r(r) / tt(S), where 
7r (^) = Sres 7r ( r )- Then, Kac's formula asserts that the mean number of steps needed to 
return to any point of E, i.e. the mean first-return time (MFRT) to X is l/7r(S). A simple 
proof of this result and of its extension to higher-order moments, which will be used later 
on, is given in Appendix IT51 
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FIG. 1: Modifications of the original lattice: arrows denote one-way links. 



Kac's formula can be used to derive the MFPT (T) by slightly modifying the original 
lattice (see FigH]): we suppress all the original links starting from the target site T, and add 
a new one-way link from T to the starting point S, whereas all other links are unchanged. 
In this new lattice, any trajectory starting from T goes to S at its first step, so that the 



MFRT to T is just the MFPT from S to T in the former lattice, plus one. 

An exact, formal expression for the MFPT can thus be derived for the most general finite 
graph. Consider iV points at positions r 1; . . . , in an arbitrary space. The transition rates 
from point j to point i are denoted Wij. If we assume that one transition takes place during 
each time unit we have: 

XX = 1 (!) 

i 

Let iy be the position of the target site, be the position of the starting site, and 7r(r) be 
the stationary probability of the modified lattice. We write 7r(ry) = J. According to Kac's 
formula, the MFRT to T on the modified graph is 1/ J, so that the MFPT from S to T in 
the original graph is (T) = 1/J — 1. All we need to find is the stationary probability n. It 
satisfies the following equation: 

7r(rj) = X w ij n ( r j) + J$iS ~ JwiT (2) 
j 

where 5 is the Kronecker symbol. To solve this equation, we define the auxiliary function 
7r', such that 7r'(rj) = vr(rj) — JtfjT- It satisfies: 

Tr'(ri) = yXjVfri) + J5 iS - J5 iT (3) 
j 

so that tt ; has the following expression: 

Tr'fo) = (1 - J)7r ( rj ) + JH(r t \r s ) - JH{Ti\r T ), (4) 
where tiq is the stationary probability of the original lattice, and H is the discrete pseudo- 

n 

Green function [23(, which satisfies the two following equations: 

HfalTj) = y^WjkHfalTj) + 5ij - — (5) 

k 

X^N ri ) = F, (6) 

i 

where H is independent of j. Moreover, if Wij is symmetric, which will be the case in all 
the practical cases considered, H will also be symmetric in its arguments. The pseudo- 
Green function can be seen as a generalization of the usual infinite-space Green function to 
a bounded domain. Indeed, Eq. <^ without the —1/N term corresponds to the definition 
of the infinite-space Green function, which would not have any solution for a finite domain 



with reflecting boundary conditions: it is necessary in this case to compensate the source 
term Sij, and the simplest way to do so is to add the —1/N term. The properties of this 
function are further discussed in Appendix [F] We can thus see that the solution (jlj) satisfies 
equation Q, and ensures that ir is normalized. 

The condition 7r'(ry) = allows us to compute J and to deduce the following exact 
expression: 

(T) = ^—[H(r T \r T )-H(r T \r s )} (7) 

If Wij is symmetric, and we will consider that this is the case in the rest of the paper, we 
simply have tiq = 1/N, and we get the simpler formula: 

(T) = N[H(t t \t t ) - H(t t \t s )] (8) 

This result may be obtained by an alternative and complementary approach. We consider 
that in the domain there is a constant flux J of particles per time unit entering the domain 
at the source point S. The particles are absorbed when they reach the target, and, since all 
particles are eventually absorbed, we have an outcoming flux J at the target. The average 
number of particles in the domain satisfies: M = J(T), which will allow the determination 
of (T). Indeed, the average density of particles p(r) satisfies the following equation: 

p(Ti) = ^ W ijP( T j) + Jhs - J$iT- (9) 

3 

The three terms of the equation correspond respectively to the diffusion of particles, the 
incoming flux in S, and the outgoing flux in T. This is exactly the same equation as Eq.(j3]), 
with the same condition p(i*r) = 0, and thus admits a similar solution, with the difference 
that the total number of particles in the domain is not fixed a priori. The solution is thus: 

p{n) = p + JH(ri\r s ) - JH( ri \r T ) (10) 

which gives, with the condition p(iy) = and the relation J(T) = M = Npo, the same 
result as before for the mean first-passage time, namely Eq.fjH]). This formula is equivalent 



to the one given in 16|, but is expressed in terms of pseudo-Green functions. One advantage 
of the present method is that it may be easily extended to more complex situations, as it 
will be shown. Another advantage is that, although the pseudo-Green function H is not 
known in general, it is well suited for approximations when the graph is a bounded regular 
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lattice. The simplest one in this case is to approximate the pseudo-Green function by its 
infinite-space limit, the "usual" Green function: H(r\r r ) ~ Gq(t — r'), which satisfies: 

G (r) = - V G (r')+5 0r . (11) 
a *■ — ' 

r'eN(r) 

where N(r) is the ensemble of neighbours of r, and a the coordination number of the lattice. 

n 

The value of Gq(0) and the asymptotic behaviour of Go are well-known [24j]. For instance, 
for the 3D cubic lattice, we have: Gq(0) = 1.516386... and Gq(t) ~ 3/(27rr) for r large. For 
the 2D square lattice, we have Gq(0) — Go(r) ~ (2/7r) ln(r) + (3/n) In 2 + 27/7?, where 7 is 
the Euler gamma constant, and (3/7r)ln2 + 27/^ = 1.029374... . These estimations of Go 
are used for all the practical applications in the following. In some cases (especially in three 
dimensions when the target is far from any boundary), approximating H by Go will give 
accurate results (see Fig. EJ). The small correction is due to boundary effects, which are 
further discussed in Section IIVI In other cases it will only give an order of magnitude. In 
the case of a rectangular or parallepipedic domain an exact expression of H is known 26], 
and the FPT from any point to any other point in the domain can be computed exactly. 
This exact result and simple approximations, which can be used in other cases, are given in 
Appendix [A] 



B. Application: absorbing opening in a reflecting boundary 

Another situation that may arise and can easily be dealt with is the case of an absorbing 
opening in a (locally) flat reflecting boundary of a bounded domain: we are interested in the 
mean time a particle takes to exit from the domain, if it may only exit by this opening, (see 
Fig. [3]). We only consider regular lattices of dimension d = 2 or 3. We can define a target 
site, just behind the flat boundary. The problem here is that the pseudo-Green function for 
the domain plus the target site is difficult to compute, whereas the pseudo-Green function 
near a flat boundary can be easily evaluated, and is even known exactly if the domain is 
rectangular or parallepipedic. To solve this problem, we will call the site next to the target 
the approach site A. We indeed have to go through this approach site in order to reach the 
target site. We will call (T)st the average time to reach the target site, starting from the 
source; (T) S a the average time to reach the approach site, still starting from the source; 
(T)^ the average time to return to the approach site, assuming the random walk does not 
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FIG. 2: (color online) 3D - Influence of the distance between the source and the target. Red 
crosses: simulations; blue dashed line: evaluation of the FPT with H = Gq. The domain is a cube 
of side 41, the target being in the middle of it. All the simulation points correspond to different 
positions of the source. 
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FIG. 3: Opening in a flat reflecting boundary 

go to the target site after exiting the approach site; {T)at the average time to reach the 
target site, starting from the approach site. We have the following equations: 



1ST 



ISA 



+ (T) 



AT 



(12) 



S 



since the random walk has to go through the approach site, and 

(Vat = ((T) AA + (T) AT ) + ± (13) 

Once the random walker is at the approach site, it may either go directly to the target site 
(probability ^, where d is the dimension of the lattice) or go another way, in which case it 
has to go back to the approach site before finding the target site. Thus, 

(T) AT = (2d-l)(T) AA + l (14) 

To compute (T)^, we have to remember that, if the boundary was fully reflecting, we would 
have the average return time: it is given by Kac's formula, and is N. We then have, with 
arguments similar to Eq. [13) 

N= 2d-1 1 

2d y 1 2d v ; 

We then have: 

(2d — 1)(T) AA = 2dN — 1 (16) 

and thus: 

(TV = 2dN (17) 

As for the average time needed to reach the approach site, starting from the starting site, it 
is exactly the same as in the case where the boundary is totally reflecting: 

(T) SA = N[H(r A \r A ) - H(r A \r s )}, (18) 

and finally: 

(T)st = N[2d + H(r A \r A ) - H(r A \r s )], (19) 

To evaluate H, we have to take into account the effect of the boundary. Since the 
boundary is flat, the simplest way to check the boundary condition is to write H(r\r') ~ 
Go(r — r') + Go(r — s(r')), where s(r) is the point symmetrical to r with respect to the bound- 
ary. We will use this approximation in the following, (cf. Appendix |A] for a discussion of 
this approximation) We note Gq(1) = Gq(0) — 1 the Green function for the sites surrounding 
the origin, and notice that T is symmetrical to A with respect to the boundary. The mean 
exit time is then: 

(T)sr ~ N[2d + G (0) + G (l) - G (r s - r A ) - G (r s - r T )] (20) 
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C. Higher-order moments 



Moreover, we are able to evaluate the higher-order moments and distribution of the FPT 
in the 3D case, provided the domain is not too elongated, i.e. the typical distance between 
a point and a boundary is N 1 ^ 3 . The computation of the moments is detailed in Appendix 
IB II However, we cannot compute the higher-order moments and distribution of the FPT 
in two dimensions, or with a too elongated 3D domain. The computational reasons behind 
this are explained in Appendix IB 1\ but we will also explain it later from a physical point of 
view. We obtain the following result for higher-order moments: 



(T n ). = n\N n [{H{v T \v T ) - H(r T \ ri )) (H{r T \r T ) - H)^ 1 + 0(nN- 2 l 3 )\ , (21) 

where H is defined by Eq. (jSJ) 

To check these results, we computed the moments with a numerical simulation (cf. Ap- 
pendix [E] for the simulation method), and found (see FigH]) a good agreement with the the- 
oretical estimation (12 II) . where H is approximated by Go, and H is approximated by is its 
value for a spherical domain, computed in the continuous limit, H = (18/5)(3/(47r)) 2//3 iV -1 / 3 
(cf. Eq. (1A10I) for the computation). 

The study of the distribution in the limit of large N will enable us to go even further. 
Indeed, if we neglect the corrections in nN~ 2 / 3 in Eq. (|2ip . the moments of T/N are those 
of the following probability density p: 

*> - ret) - y ) - (-wj) + is^i '«> < 22) 

The large- N limit of this probability density is rigorous (since the corrections to the moments 
vanish). In this limit, if (r^rr) tends to Gq(0); H tends to 0. Thus, the probability density 
of T/N tends to the following probability density, the relative position of % and T being 
fixed: 

m = ( gW) ) exp {-gm) + ^<fir m (23) 

These results have been confronted to numerical simulations (FigJSJ) . We computed the 
exact distribution for several domain sizes, and may notice that the curve divides in two at 
short times. This is due to the fact that, at short times, the parity of the step is important: 
as long as the walk does not touch the boundary, the distance between the starting point 
and the walker has the same parity as the time elapsed. The time needed for the two curves 
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40 

FIG. 4: (color online) 3D - Relative difference between the simulations and the theoretical prevision 
(|2ip . The domain is a cube of side 51 centered on the target at (0,0,0) and the source at (2,2,1). 
The order of magnitude of the relative difference is indeed nN~ 2 /^ . 

to collapse into one shows very well the time needed to erase the memory of the starting 
position. The curves before this time correspond to the Dirac part of the probability density 
f l22|) ; however, we can see that, once the two curves have collapsed, the resulting curves 
fit very well the theoretical prediction (|22l) . which is indeed more accurate than the limit 
probability density ( 1231) . 

To analyse the physical meaning of this result, we may first notice that, if the probability 
density (1221) is averaged on the starting point, the Dirac part of the density vanishes, and 
we simply have an exponential distribution of the first-passage time. This property sheds a 
new light on the quasi-chemical approximation [3], which assumes that if a particle starts 
randomly in a volume, and may only exit through a small hole, it has a constant probability 
to exit at each time step. This approximation leads to an exponential distribution of the 
exit times. If we consider that the target site is the exit point for a particle, then the exit 
time is exactly the FPT. Thus, we have an evaluation of the accuracy of the quasi-chemical 
approximation (or at least of its moments) in this case. 

The interpretation of the probability density ff22l) is the following: the first part of the 
density, which decays exponentially, corresponds to the decay of the probability distribution 
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FIG. 5: (color online) Simulation of the probability distribution of the FPT, rescaled as a proba- 
bility density, with different domain sizes. In both cases, the target is in the middle of a cube at 
position (0,0,0) and the source is in (1,2,2). We plot the estimated density (|22p vs. numerical sim- 
ulations for different domain sizes. The dark blue (simulation) and cyan (estimated density) curves 
correspond to a cube of side 11 (N = 1331); The red (simulation) and orange (estimated density) 
correspond to a cube of side 21 (N = 9261) ; The green dashed curve corresponds to the high-iV 
limit (|23p of the probability density. For both domain sizes, the simulated distribution splits into 
two parts for short times and cannot be distinguished from the theoretical density afterwards. The 
labels a and b correspond to domain sizes of 11 and 21, and the high-iV limit is labeled c. 

of the FPT if the particle starts randomly in the set of points. The second part corresponds 
to a particle reaching the target in a time negligible with respect to N. Here we must 
remember that a free 3D walk is transient: the particle may never reach the target in 
infinite space. Thus, one can interpret the Dirac term as the probability to reach the target 
without touching the boundaries. For N large enough it is equivalent to the probability to 
reach the target at all in infinite space. And, for this kind of trajectory, the probability 
distribution of the FPT does not depend of N, and, thus, the probability density of T/N 
will tend to 5(0) for large N. On the other hand, if the particle does reach the boundary 
(it happens after a typical time iV 2//3 , since the boundaries are at a typical distance N 1 ^ 3 , 
and the typical time needed to cross a distance r is r 2 ), its position will become random in 
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FIG. 6: (color online) Simulation of the distribution of the FPT in the 2D case, rescaled as a 
probability density. The domains are squares of different size; The target is in the middle of the 
squares at position (0,0) and the source is in (2,3). The side of the squares are: 21 (red curve [a]), 
41 (orange curve [b]), 81 (yellow curve [c]) and 161 (green curve [d]). The semi-logarithmic scale 
shows the long-time exponential decay. The splitting of the curves for short times is due to parity 
effects. 

a time negligible with respect to N, and, thus, the probability density of T/N will be the 
same as if the particle started in a random position in this latter case. This argument fails 
for an elongated domain, which can be seen as a physical reason why we are not able to 
compute the FPT distribution in this case. We can check that the probability to reach the 
origin for a random walk in infinite space is indeed q°J^ 24. 25]. 

Here we can see an important physical difference between the 2D and 3D cases: in two 
dimensions the random walk is recurrent. We can thus conclude that the large- iV limit 
probability density of T/N will be a simple delta function, since, in the limit of infinite 
space, the particle almost certainly reaches the target in a finite time, even if the MFPT 
is infinite! However, the probability distribution for finite N will be much more difficult 
to compute: Indeed, there will also be two regimes, of low T/N, when the particles have 
not touched the boundary and thus the distribution is the same as in free space; and the 
regime of high T/N, where the distribution decays exponentially (since the system has lost 
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the memory of its starting point). The transition between the two regimes happens at a 
finite T/N (since the time needed to reach the boundaries is of order N). Thus, the low 
T/N regime will have a much stronger influence on the values of the moments than on the 
3D case, which may explain why the computation of the moments and distribution is much 
more delicate in this case. We can see in Fig. [H] typical probability distributions for different 
domain sizes. One can very well see that the transition between the two regimes takes place 
at a finite T/N no matter the size of the domain, and that the long-time regime indeed 
corresponds to an exponential decay. 



D. Case of two targets 

We can now assume that the lattice contains not one but two target points T\ and T 2 . 
The problems that may arise in this case are the mean time needed to reach one of the two 
targets, which we will call mean absorption time and note (T), and the splitting probabilities, 
i.e. the probabilities Pi to reach T\ before T 2 and P 2 to reach T 2 before T\. This model 
corresponds to the case of a diffusing particle which may be absorbed either by the target 
Ti or the target T 2 . We can also, even if it will be less straightforward, study the conditional 
mean absorption time, i.e. the mean absorption time (Ti) (resp. (T 2 )), for particles which 
are absorbed by the target T\ (resp. T 2 ). This is relevant in many chemical applications 
[2I, and may be useful in^biology to determine to which extent cellular variability may be 



controlled by diffusion [271 ]. 



To compute these quantities, it is more convenient to use the alternative approach pre- 
sented on page [6j we consider a constant incoming flux of particles J, and we have an 
average outcoming flux of particles J\ in T 1; and J 2 in T 2 . Since all particles are eventually 
absorbed, J\ + J 2 = J. The probability to reach the target i is then Pj = Ji/J. The total 
number of particles M in the domain satisfies M = J{T), and the mean density of particles 
satisfies the following equation: 

p(r;) = ^ Wjjpjrj) + JS iS - JiS iTl - J 2 6 m (24) 
j 

We then get: 

p(n) = Po + JH(n\r s ) - J!H(ri\r Tl ) - J 2 H(vi\r T2 ), (25) 
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then, writing p^tJ = p(fT 2 ) = 0, we get the following set of equations: 



po + JH ls - JPiHo, 
< p + JH 2s - JP 2 H 02 



JP 2 H 



JP l H 



12 — 



12 — 











(26) 



P1 + P2 



1 



V 



where H 12 = H(r Tl \r T2 ) and, for i — 1 or 2, H is = H(r Ti \r s ), H 0i = H(r Ti \r Ti ). From this 
set of equation we can deduce Pi, P 2 and po = J(T)/N. We thus get exact expressions for 
the mean absorption time and the splitting probabilities, respectively: 



This result can be extended if necessary to more than two targets; if there are n targets, 
we have n+1 unknown variables (po and the n probabilities P&), with n+1 equations, namely 
^2 Pk = 1 and the n equations p(r*r fc ) = 0, which is enough to determine all the unknown 
variables. However, this may quickly become computationally expensive for a large number 
of targets. 

We compared the two-target results to simulations (Fig. [7]). Note that if we use the exact 
value for H, which we can compute for a cube (cf. Appendix [AJ , it is indeed impossible to 
see a difference between the theoretical predictions and the simulations. 

It is interesting to underline an important qualitative difference between the 2D and 3D 
cases. In 3D, the furthest target always has a significant probability to be reached first, 
since the most important terms in the probabilities Pj are H m and H 02 . In 2D, if a target 
is much closer from the source than the other, it will almost certainly be reached first, since 
H(ri\rj) scales as In |r^ — r j | . Actually, the probability for the furthest target to be reached 
first decreases logarithmically. These properties are related to the transient character of 
the infinite 3D walk, and the recurrent character of the 2D walk: indeed, an infinite 2D 
walk explores all the sites of the lattice, whereas an infinite 3D walk does not; we may thus 
consider that the 2D walk will explore most of the sites surrounding the source before going 
much further, whereas the 3D walk will not, which qualitatively explains the difference of 
behaviour. 

We can also determine the conditional absorption times (Ti) and (T2). For this, we 
will compute A4, the average number of particles in the domain which will eventually be 



(T) = N 



(Hpi — Hls){H<i2 — Pl2s) — (Hl2 — Pl2s){Hl2 ~ His) 

Hqi + H 02 — 2H 12 



(27) 




(28) 
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FIG. 7: (color online) 3D: two-target simulations. Simulations (red crosses) vs. theory with the 
approximation H = Gq (solid line). One target is fixed in (-5,0,0); the source is fixed in (5,0,0); 
the other target is at (x,3,0). The domain is a cube of side 51, the middle is the point (0,0,0). 

absorbed by TV We have A4 = Jfc(Tfc), which will allow us to compute (T^). To compute 
Afc, we can simply notice that the density of particles that will eventually be absorbed by 
Tk at the point i is simply p(rj)Pfc(rj), where Pk{^i) is the probability to be absorbed by Tk 
if the walk starts from i. We thus have: 

A4 = ^p(ri)P fe (rO (29) 
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This equation is exact but may prove quite difficult to compute, especially in two dimensions 
if H is not known exactly. However, in three dimensions, we may use the same kind of 
approximations as for the computation of the high-order moments of the FPT (with the same 
limitations, i.e. the 3D domain should not be too elongated) to estimate the conditional 
probabilities. If we note H^s = if(rj|r s ) and = H(ri\rx k ), we have: 



\r (Ha — H i2 + H 02 — Hi 2 )(po + J H iS — J\Hn — J 2 Hi2) , , 

• Nl ~ 1^ u T tt KTT l dU J 



We use the properties J2i H ( r i\ r j) = NH ( cf - E Q- ©) and J2i H { r i\ r j) H { r i\ r k) = 0(N 1 ^) 
(cf. Eq. IET21) to write: 

m = N (_MB^ki±m^i (31 ) 

And we can conclude: 

(Tl) = i- g " g -^ + y <T> (32, 

"l -"01 + n 02 — £n 12 

The expression for (T 2 ) is of course equivalent. This expression is not exact, but is very 
accurate: the relative difference between the numerical simulations (see FigJH]) and the 
expression ( 1321) is of about 0.01%, for a domain of size = 51 3 . 

Finally, we have a wide range of quantities which can be computed exactly, or with a 
very good accuracy, provided we know the pseudo-Green function H. Unfortunately, there 
are only a few cases in which it can be computed exactly. Otherwise, we will have to 
use approximations, which, of course, give less accurate results. Both exact results and 
approximations are detailed in Appendix lAl 



III. BROWNIAN MOTION ON CONTINUOUS MEDIA 

We may consider a similar problem in a continuous medium (see Fig. [H]): if we have a 
Brownian motion whose diffusion coefficient is D, how much time does it take to reach a 
target? A difference with the discrete case is that the target has a finite size a which is 
an important parameter of the problem. We will consider a spherical target T, of radius 
a, centered in r^. The Brownian motion starts from the starting point 5* (its position is 
denoted by rg). It is restricted to a domain T> of volume V (for 2D domains we will call 
the area A) , and we note D* the domain deprived of the target. We will derive the same 
quantities as in the discrete case, but the results are this time only approximate; we can thus 
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FIG. 8: 3D: two-target simulations. The conditions are identical to those of FigJTJ we show the 
conditional absorption times (Ti) (resp. (T2)). The blue crosses (resp. red plusses) show the 
results of the numerical simulations, the cyan [a] (resp. orange [b]) dashed line shows the theoret- 
ical expression (f3"2j) with H = Go, the green [c] (resp. brown [d]) solid line shows the theoretical 
expression with the exact value of H ()A3|) . 



add some refinements to the method, in order to increase the accuracy. These refinements 
are given in Appendix [Cj and are used in practical computations of the MFPT in Appendix 
[A] when the target is close to a boundary. It should be emphasised that, in the cases where 
the pseudo-Green function is known, such as the case of a spherical domain, the method 
gives accurate explicit expressions for all the MFPT and the other quantities studied here. 




FIG. 9: (color online) Continuous problem 
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A. Mean first passage time 



The mean first passage time (MFPT) (T(r s )) at the target satisfies the following equations 

DA(T(r s )} = -1 if r s G V* (33) 

(T(r s )> = if r s G S abs (34) 

0„(T(r s )> = if r s G S refl (35) 

where S abs (resp. S rc fl) stands for the surface of the absorbing target sphere (resp. the 
reflecting confining surface) and d n denotes the normal derivative. The boundaries have to 
be regular enough (twice continuously different iable is sufficient, but not necessary) for these 
definitions to make sense. 

To solve this problem, we introduce the following Green function G(r|r') defined by 

- AG(r|r') = 5{r - r') if r G V* (36) 

G(r|r') = if r G S abs (37) 

<9 n G(r|r') = if r G S rcfl (38) 

Note that this Green function may also be seen as the stationary density of particles if there 
is an unit incoming flux of particles in r', and the diffusion coefficient is set to 1. It should 
not be confused with the free Green function Go, and is rather the continuous equivalent of 
the average density of particles p defined in Eq. with J = 1. It depends implicitly on 
the target position through Eq. (T3T|) . 
Using Green's formula, 

/ ((T(r))AG(r|r')-G(r|r')A(T(r)))^r= j «T(r)> <9 n G(r|r') - G{v\v')d n (T(r))) d^r, 

(39) 

we easily find that the MFPT is given by 

(T(r 5 )> = 1 G(v\v s )d d v (40) 

To approximate G(r|r5) we can use a direct transposition to the continuous case of Eq. 
TDD : 

G(r|r s ) ~ p (r s ) + H(t\t s ) - H(t\t t ) (41) 
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where po is defined by G(t\ts) — if r G S a bs and H(r\r') is the pseudo-Green function [231 ]. 
which satisfies: 

- A#(r|r') = 8(r - r') - if r <E V (42) 
<9 n #(r|r') = if r G £ rcfl (43) 
i?(r|r') = #(r'|r) (44) 
H(r'\r)d d r' = VH, (45) 



v 

H being independent of r. This latter equation can be easily deduced from the three previous 
ones. 

The choice ( 14T|) of G(r|r') is the simplest one which satisfies formally Eqs ( l36l) and 
(1381) . However, (1371) can only be approximately satisfied. To take into account this latter 
equation, we will approximate, on the target sphere, if (r|rs) by H(yt\ys) and H(y\yt) by 
Go(r — r*r) + H* (iy|iy) , where Go is the well-known free Green function( (2tt)^ 1 ln(r) in 2D, 
l/(47rr) in 3D), and H* is defined by: 

H*(r\r')=H(r\r')-G (r-r'). (46) 

Note that if*(r|iY) has no singularity in r^. Thus on the surface of the target sphere we 
have: 

p (r s ) + H(r T \r s ) - G (a) - H*(r T \r T ) = 0, (47) 
where Go (a) is the value of Go(r) when |r| = a. We can now compute 

(T(r 5 )> = ^J (Po(r s ) + H{r\r s ) - H{r\r T )) d d v (48) 

Since the target is small compared to the domain, the integral over T>* is almost equal to 
the integral over D, the relative order of magnitude of the correction being a 3 /V in 3D and 
a 2 1 A in 2D. Using the property f|45|) . we can then compute the integral, and find the result: 

<T(r 5 )> = = I (G (a) + H-frM - H(r T \r s) ) + O ( ^M) (49 ) 



D D y *i °» ■ y D 

This equation is very close to ([8]), with the correspondence H(r\r) — > G (a) + H*(r\r), 
but one should pay attention to the fact that this is only an approximation! One may 
expect deviations from this expression when the variations of iJ(r|rs') or H*(t\tt) will not 
be negligible over the target sphere; it corresponds to the cases when the target is either 
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FIG. 10: (color online) Brownian motion on a 2D disk of radius 25 centered on (0,0); the source is 
in (0,1) and the target of radius 1 is on (x,0). Red crosses: simulations; black solid line: estimation 
l9j) with the exact function H for a sphere given by Eq. (|A6|) 



near the source or near a boundary. However, if we compare the expression obtained with 

simulations (see Fig JTUi) when the target is near the source, we see no such deviation; this is 

justified in Appendix O On the other hand, there is indeed a deviation near the boundaries. 

This deviation scales as a/d in two dimensions, or a/d 2 in three dimensions, where d is the 

distance between the target and the boundary. It is possible to compute a correction, which 

is explicited in Appendix EH and used in practical applications in Appendix [A] 

The exact value of H is known analytically for disks and spheres [23], we will detail this 

in Appendix [A] This is why we will test the expressions we obtain in such geometries. If 

no exact expression is known, the simplest approximation of H is simply H = Gq. More 

accurate approximations are also discussed in Appendix |A] We give the estimations of 

(T(rs)) with the basic approximation, to show the order of magnitude: 

V (\ 1 
4^D \a~R 



(T(r s )> 



(3D) 



(50) 



(T(r s )) 



— In- (2D) 



(51) 



R being the source-target distance. This already improves the (exact) asymptotic results of 
Pinsky j^j, which only give the leading term in a. 
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B. Higher-order moments 

The higher-order moments and density of the FPT in the three-dimensional case can also 
be computed. The computation is detailed in Appendix IB 21 the results are quite similar to 
the results obtained in the discrete case, and the physical interpretation is essentially the 
same. The results obtained are the following: 



(T"(r 5 )> = [(G (a) + H*(v T \r T ) - H(r T \r s )) (G (a) + H*(v T \r T ) - H) + O («r 2 / 3 a 2 "") 

(52) 

We may also deduce from this information about the probability density of the absorption 
time p(t): If we drop the term O [nV~ 2 ^a 2 ~ n ^j , we have: 

Pit) = D G (a) + H*(r T \r T )-H(r T \r s ) ( -Dt \ 

V (G (a) + H*(r T \r T ) - S) 2 1 \V (G (a) + H*(r T \r T ) - H) J 
■ H(v T \r s )-H 
G (a) + H*(r T \r T ) - H 1 ' 
In the limit a — > 0, with the position of fixed, the H terms are constant since they 
only depend on the shape of the domain, and Gq(o) tends towards infinity. The probability 
density then simply becomes exponential: 

AnaD ( 4naDt\ 



p(t) = -^exp( — ) (54) 



In the limit a — > 0, with the quantity R/a fixed, H(ts\tt) ~ Gq(R), and the probability 
density becomes: 

, , An Da ( a\ ( AitaDtX a „. , . 

P(t) = — (i- 5 )«p(-— J+ 5 *M < 55 > 

We did not test these results with a numerical simulation, since the continuous simulation 
method (see Appendix [Ej) is not adapted to the computation of the FPT density, and would 
require a large computation time to give accurate results. Furthermore, the approximations 
made (cf. Appendix [Bj are the same as on the discrete case, and the discrete results have 
been successfully compared to an exact numerical simulation (cf. Fig. [5]). 



C. Case of two targets 

For the case of two targets, we will compute the same quantities as in the discrete case; 
however, we may notice that the radius a x and a 2 of the two targets may differ, which adds 
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another parameter to the problem. With two targets, we will use the same Green function 
as before, but £ a b s = £1 + £ 2 will be the reunion of the surfaces of the two absorbing 
target spheres. The mean absorption time (T(r s )) satisfies the equation (HUj) ; the splitting 
probability Pi(rs) satisfies the following equations 

APx(r) = (56) 

Pi(r) = 1 if r G Ei (57) 

Pi(r) = if r G £ 2 (58) 

d n P x (v) = if r G S refl (59) 

Using Green's formula, we get: 

P 1 (r s ) = - f d n G(r\r s )dr, (60) 



The expression for P 2 is of course similar. Note that the normal derivative is oriented towards 
the inside of the target. A simple approximation of G, equivalent to the discrete Eq. f[2"5"j) 
is: 

G(r|r s ) = Po(rs) + H(r\r s ) - P 1 (r s )H(r\r Tl ) - P 2 (r s )Pf (r|r T2 ). (61) 

This expression satisfies Eq. (j3"5"|) .( }3"8"|) and flBDl) . and p , Pi and P 2 are set in order to satisfy 
Eq. (l3"7j) approximately. We use the same approximations as in the one-target case, which 
gives the following set of equations: 

p (r s ) + H u - PiPqi - ^2^12 = 

Po(r s ) + H 2s - P 2 H m - P x H l2 = (62) 

Pi + P 2 = 1 

where H 12 = P(r Tl |r T J and, for i = 1 or 2, H is = H(t t .\t s ), H 0i = G (ai) + H*(r Ti \r Ti )- 
These equations are exactly identical to the discrete equations, only the meaning of the H 0i 
changes. We thus can deduce, using the same relation between p and (T) as in Eq. 



(T(r 5 )) — — ^ 01 - ^ ls ^^ 02 ~ ^ 2s ^ ~ ~ ^)(Hi2 - Hi s ) ^gg^ 
D Hqi + H 02 — 2H 12 

_ H ls +H 2-H 2 s-Hv2 

H 01 +H 2-2H 12 ^g^-j 

_ H2s+Hgi—Hi s —Hi2 

Hoi+H 02 -2Hi2 
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FIG. 11: (color online) Brownian motion on a 2D disk of radius 25 centered on (0,0); D=l; the 
source is in (-5,2) and the two targets of radius 1 are on (5,2) (Ti) and (x,0) (T2). Red crosses: 
simulations; black solid line: estimations (|63[) and (|64p with the exact function H for a disk (|A6p . 

We show in Fig. (TT] and [12] the results of the numerical simulations. We can see that they 
are accurate, with a small correction (the relative correction scales as (a/d) \n(d/a) in 2D or 
a 2 /d 2 in 3D, d being the distance between the two targets) when the two targets are close 
(an explicit correction is given in Appendix [C]), and when one target is near a boundary 
(exactly as in the one-target case). 
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FIG. 12: (color online) Brownian motion on a 3D sphere of radius 25 centered on (0,0,0); D=l; the 
source is in (-5,2,0) and the two targets are on (5,2,0) (Ti, of radius 0.5) and (x,0,0) (T2, of radius 
1.5). Red crosses: simulations; black solid line: estimations (j63j) . (|64[) . with the exact function H 
for a sphere (|A7p . 

The curves themselves deserve a few qualitative remarks. Unsurprisingly, the splitting 
probability P2 is maximal when T 2 is the closest to the source. When the two targets have 
different sizes, an interesting phenomenon appears (Fig. [T2|) : the probability to hit the 
largest target (T2) has a second maximum when it is close to the other target. One can 
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understand this by a scaling argument. If the two targets are far away, Pi will be about 
Oi/(ai + 0.2)- If the two targets touch one another, and a% <C ci2, then the target T\ covers 
a surface iraf of the target T 2 . It can thus be expected that the probability Pi will scale as 
a\/a\, and thus be much lower than if the two targets were far away These arguments are 
for the 3D case, but the qualitative behavior would be the same in the 2D case. However, 
the behavior of the splitting probabilities when one target is much further than the other 
from the source will be different in 2D and 3D, for the same reasons as in the discrete case. 
In the figures the domain is not large enough to make the difference obvious. The mean 
absorption time has a similar qualitative behavior in both cases: an unsurprising minimum 
when the moving target is close to the source, maxima when the moving target is near a 
boundary, due to boundary effects, and a maximum when the two targets are close, which 
deserves a few more comments. This could indeed be predicted directly from Eq. fl63|) . but, 
physically, this comes from the fact that, if the two targets are close, a particle undergoing 
a Brownian motion, which reaches one target, often would have reached the other shortly 
afterwards in a single target situation. Thus, the mean time gained, compared to the single 
target situation, will be much lower when the two targets are close. To analyse the values 
themselves, one should keep in mind that the times are normalized by V/D\ the order of 
magnitude of the normalized times will then be Go(a) — Gq(R), which explains the values 
around 0.05 obtained in the 3D case. 

As for the conditional FPTs (Ti(rg)) and (T^r^)}, we have the following relations [l|: 

DA(Pi(r)(Ti(r))) = -P x (r) if r G V* (65) 

P 1 (r)(T 1 (r)) = 0ifrGS abs (66) 

^(P 1 (r)<T 1 (r))) = 0if reS teJl (67) 

and of course the equivalent relations for (T 2 (r)). We use as usual Green's formula, and 
obtain: 

P 1 (r s )(T 1 (r s ))= / G{v\v s )P 1 {v)dv (68) 
Jv 

This equation is very similar to the discrete Eq. (129]) and the following calculations for the 
3D case are exactly identical, and give: 

<Tl(rs)> = P,(r s) H m+ H m -2H 12 <T(rs)> (69) 
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FIG. 13: (color online) Two-target simulations. The conditions are identical to those of Fig. [12\ 
we show the conditional absorption times (Ti) (resp. (T2)). The blue Xs (resp. red +s) show 
the results of numerical simulations; The greenfa] (resp. brown[b]) solid line shows the theoretical 
estimation (|69p with the exact function H for a sphere (|A7|) . 

We show in Fig {13] the result of numerical simulations. The noise is more important than 
in other simulations, especially for (Ti). This is due to the fact that the probability Pi is 
often small, which reduces the number of processes on which the time is averaged, and thus 
increases the noise. 

We thus are able to compute first-passage times, splitting probabilities and absorption 
times with a good accuracy (especially with the improvements given in Appendix [C]), pro- 
vided we know the pseudo-Green function H. The computation of H is discussed extensively 
in Appendix [A] and more briefly in the following. 

IV. DISCUSSION 

The computation of the pseudo-Green function can be a difficult problem. Indeed, there 
are a few cases when it can be computed exactly (see Appendix IA II) , namely in the discrete 
case for a rectangular/parallepipedic domain or for periodic boundary conditions, and in the 
continuous case when the domain is a disk, a sphere, or the surface of a sphere. Otherwise, we 
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have to use an approximation, the simplest ones being presented and discussed in Appendix 



In the following we present a synthetic and qualitative description of the important pa- 
rameters which have to be taken into account when it comes to computing the mean-first 
passage time. 

The first and most important parameter is the size of the domain. Indeed, the MFPT 
is proportional to the size of the domain, both in two and three dimensions. The second 
essential parameter is the size of the target for the continuous case: once we have these two 
parameters we already have a rough order of magnitude of the MFPT. The third important 
parameter is the distance between the source and the target. In three dimensions this 
parameter is important as long as it is of the same order of magnitude as the target size; 
its influence is inversely proportional to the source-target distance. In two dimensions this 
parameter will be important at any distance, since the MFPT depends logarithmically of this 
distance. Once these parameter have been taken into account we have a good approximation 
of the MFPT (Eq. (1511) and (ED}) if both source and target are far from any boundary. To 
see what far means in this case, a good criterion is that any correction involving a boundary 
(see below) is negligible. Otherwise we have an order of magnitude, and to proceed further 
we will have to take into account the precise position of the boundaries. 

The qualitative effect of the boundary is to increase the MFPT when the target is near 
a boundary, and to decrease it when the source is near a boundary (it can be seen in the 
following equations). The first effect is much more important than the second: in three 
dimensions, with a flat boundary, a basic approximation gives: 



where s(r) denotes the point symmetrical to r with respect to the boundary. One can see 
that the influence of the boundary is inversely proportional to the distance between the 
target and the boundary. This is also true if the source is near a boundary, which is why 
the most important parameter is indeed the position of the target. One may note, however, 
that if the target or the source lies in a corner, these effects are amplified. 

In two dimensions the influence of the position of the boundary is more important, and the 
position of the source is a relevant parameter: a basic approximation with a flat boundary 




(70) 
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gives: 

(T) = — f In |rg ~ rT| + In ^4^4 1 (71) 
2vr V a |s(r T ) -r T |/ v ; 

If the target is much closer to the boundary than the source the effect can be to double the 
MFPT; on the other hand, if the source only is near a boundary, the related correction is 
bounded. The corners also have an amplifying effect in two dimensions. 

The quantitative estimates thus obtained are generally more accurate in three dimensions 
than in two dimensions, due to the fact that the effect of the boundaries on the pseudo-Green 
function is essentially local in three dimensions. In two dimensions there is still room for 
improvement, but an extensive discussion would be beyond the scope of this article. 

Conclusion 

In this article we managed to compute the mean first-passage times, the splitting proba- 
bility and the full probability density of the first-passage time (in three dimensions) with a 
good accuracy for spherical or rectangular domains. For other shapes (with a regular enough 
boundary), we gave the basic tools to approximately estimate these quantities. These results 
are especially important in the analysis of diffusion-limited reactions: The first-passage time 
corresponds to the reaction time if one of the reactants is static, and the reaction rate is 
infinite. Two promising extensions of our work would be to take into account finite reaction 
rates, which would increase the relevance of our work to reaction-diffusion processes; and 
to study the same problem with anomalous diffusion, which is relevant in many physical 
situations. 
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APPENDIX A: EVALUATION OF THE PSEUDO-GREEN FUNCTION 



1. Exact formulas 



a. Periodic boundary condition and rectangular domains for a discrete pseudo-Green function 

There are two specific cases where the discrete pseudo-Green function H may be computed 
exactly: when the domain is rectangular (parallepipedic in three dimensions) or when the 
boundary conditions are periodic [26|. These results are interesting in themselves, but, 
moreover, for a domain which is almost rectangular/parallepipedic, they will give a good 
approximation for H. For periodic boundary conditions, if we consider a domain with 
X sites in the x direction, Y sites in the y direction, and Z sites in the z direction, a 
straightforward Fourier analysis gives: 

X-l Y-l Z-1 ov -r. ( 2imn(x-x') j_ 2inw(y-y') _j_ 2ipTr(z-z') 



N±^t^ n . 1- I (cos^ + cos^ + cos^) ( 1} 



m=0 n=0 p=<5( m ,„) (0 ,o) 

In two dimensions, we have a similar formula for H: 



X-l Y-l 



2irmr(x—x') . 2imr(y—y') 



H(r\v>) - 1 V T (A2) 



exp ( ± 1 + 

... (cos ^ + cos ^) 

m=0 m=S n0 2 \ X Y I 

For a parallepipedic domain we get a slightly more complicated expression, and we have 
to use semi-integer coordinates for the points: x(resp. y and z) varies between 1/2 and 
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X(resp. Y and Z) —1/2. The result is the following: 

. , 8 cos 2^ cos cos cos 2M* cos ^ cos 2^ 
ff r r' = — > > > r-r - ^ — (A3) 

m=l n=l p=l o V A r / 

6 ^^cos^cos^cos^cos^ 6 ^cos^cos^ 
+ l-I(cos^ + cos^) + N ^ 1 - cos f 

m=l n=l ^ V A r / p=l ^ 



6 cos cos cos 2|£ cos £|£ 6 ^ cos ^ cos HSK 

+]V 2-j 1 - 1 (cos + cos + iV ^ 1 - cos f 

m=l p=l ^ V A Z / n=l 

V— 1 2—1 n7rt/' p7rz' mry p7rz « X— 1 mux' mnx 



6 y^ y^ COS COS ^-g- cos y cos z i 6 ^ COS — |— cos v 



iV^^ l-|(cos^ + cosf) iV 2 — 1-cosf 

n=l p=l ■£ \ * ^ / m=l & 

In two dimensions the expression is slightly less imposing: 



4 ^ K^ 1 cos cos COS 2jp cos ^ 



iJfrlr') = — \ \ x — x 

N ^ ^ 1 - 1 (cos ^ + cos ^) 

m=l n=l ^ v A r / 

4 cos 2^ cos 2f* 4 cos ^ cos ^ 

+ iV^ l-COS^ + iV^ l-COS,? ( ^ 

m=l a n=l * 

These formulae have the advantage of being exact, which enables us to compute exactly 
all the quantities studied in this article for such geometries. However, the computation of 
H may be computationally expensive for large domains. In the continuous case, the same 
method can be applied, but H can only be expressed as an infinite series |23l|. We give the 
result for a 2D rectangle X x Y: 



4 ~ ~ COS ^ COS COS ^ COS ^ 



XY ( mn \ _i_ ( sol) ' 

m =l n=l \ X ) \ Y ) 



2 ^ CQS2^CQS2H£ 2 ^ COS^COS^ 
jy 2-/ (-miry xy 2-^ (HzTp 



m=l 



b. Disks and spheres for the continuous pseudo-Green functions 

In the continuous case there is however a case where the pseudo- Green function is known 
exactly: if the domain is a disk or sphere of radius b. We will simply give the results; the 
detailed computation can be found in [23[. In both formulas, we use the image of r', that 
we note ?', which is aligned with r and the center of the disk/sphere O, and at a distance 
r' = b 2 /r'. We note R = |r — r'|, R = \r — f'\, and /x = C0S7, 7 being the angle between r 
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FIG. 14: (color online) Schematic picture of the quantities used in the computation of H(r\r'). 

and r'. In two dimensions, the result is the following: 

1/6 6 6 r 2 + r' 2 \ 

ff < r ' r '>^( ln S + ln 7? + 1 V + ^H (A6) 

The first term corresponds to Gq, the second to the image of r', the third term is needed to 
ensure the symmetry of H, and the last term corresponds to the — 1/V term in the definition 
of the pseudo-Green function. The three-dimensional result is a bit more complicated, with 
a logarithmic term whose physical signification is unclear: 

m i a 1 ( 1 b 1, fr'R n rr'A r 2 + r ,2 \ 
H(t\t') = — - + ^ _ in — + 1 - —f + A7 
v 1 7 4tt \ R r 'i2 6 V 6 2 6 2 / 26 3 y v ' 

These results are very useful by themselves, but they will also be useful to approximate H 
near a curved boundary, as we will see in the following. The result for a sphere can also be 
used to estimate H when one uses the approximation H = Gq in non-elongated 3D domains. 
Indeed the exact result enables one to take into account the corrections to Go, which are 
negligible when the source and the target are close, but give a substantial correction to the 
value of H. To compute H, one can use Eq. (145j) . and choose for r the centre of the sphere. 
We have in this case: 

*< r 'i r = >-i(5 + S) < A8 > 
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A constant (1 — ln(2))/(47r&) has been suppressed, in order to have a final result relevant 
for the approximation H = Gq. From this expression of H it is straightforward to get an 
expression for H: 

H = Ul\ 23 V -i/3 (A9) 



H = H (JL ) AT-i/3 (A10) 



5 \4tt J 

If one wants to use this result in the discrete case, it should be noted that the continuous 
limit of the discrete model corresponds to D = l/2d and not D = 1. This diffusion coefficient 
is included in the discrete pseudo-Green function, and the discrete estimation of H is thus: 

2/3 

5 V47T 

c. Surface of spheres 

Another case where we can compute exactly H is the case of the surface of a sphere. 
Indeed in this case we have exactly: 

#(r|r') = — -ln|r-r'| (All) 

Since H is isotropic in this case it simplifies things: G (a) + H*(r T \r T ) can be replaced 
by H(a) in Eq. (1491) This gives back the result obtained by a straightforward computation 



of the FPTs in a sphere 29j. Moreover this will give good approximations of all the two- 
target quantities, which was, to our knowledge, not known until now. This result is not 
used elsewhere in the paper, but is however important due to the physical relevance of the 
diffusion on the surface of a sphere. 



2. Use of the approximations 

The next step is to study cases where no exact formula for H is known. The simplest 
approximation to H is the infinite-space Green function Gq, but this approximation in often 
unsatisfying. We thus present a few ways to improve it. Before we present them, it must be 
emphasised that, in general, all the H terms should be derived with the same approxima- 
tion: H is defined up to a constant, and this constant depends of the approximation used! 
However, for complicated expression involving H, this constraint can be relaxed: if the ex- 
pression can be decomposed into terms of the form (H(ri\r2) — Hlr^lr^)), these terms may 
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be computed with different approximations, since they do not depend on the constant up to 
which H is defined. For example, in the two-target problem, we have Pi = ■ 
We can use if necessary two approximations, one accurate around Ti, which we note 
and another accurate around T2, which we note H^ 2 \ Then, to compute P±, we use them 
the following way: 



p _ -"Is 1 -"U2 -"2s -"12 ( \-\o\ 

H {1) + H {2) - H {1) - H {2) 1 J 

-"01 ^-"02 11 12 n 12 

This trick can be especially useful if one has to deal with two targets near two different 
boundaries. 



3. Approximations 

The most basic approximation is already known: it is the approximation H = G$. Its 
physical meaning is to ignore the presence of the boundaries, as far as the pseudo-Green 
function is concerned. To improve this approximation, there are essentially two ways: the 
first is to take the boundaries into account locally, and to satisfy the boundary conditions 
at the nearest boundary, we will see how in the following. The second one is to take the 
boundaries into account globally, by taking into account the terms —1/N or —1/V in the 
definition of H. The order of magnitude of the related correction will be of about (r — r') 2 /N 
in the discrete case, or (r — r') 2 /AA in the 2D continuous case, (r — r') 2 /QV in the 3D 
continuous case. It is thus much weaker in 3D (the maximal relative correction scales as 
jy-1/3 or a /y 1 /^ than in 2D (where the maximal relative correction scales as l/ln(iV) or 
l/\n(y/a 3 )). A more detailed discussion of this kind of corrections would be technical and 
beyond the scope of this article, but the above order of magnitude can be a good evaluation 
of the accuracy of the following boundary approximation. 

This approximation takes explicitly into account a planar boundary, and ignores all the 
others. It can be used both in the continuous and in the discrete case. If we note s(r) the 
point symmetrical to r with respect to the boundary, then the local approximation: 

H(r\r') = G (r - r') + G (r - s(r')) (A13) 

satisfies the boundary conditions on the flat boundary, and is symmetric. It thus can be a 
good approximation for the pseudo-Green function. 
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FIG. 15: (color online) Discrete random walk - Influence of the position of the source; the domain is 
a 2D square of side 41 centered on (0, 0, 0), the target is at (—18, 0) and the source is at (x — 20, 0); 
blue dashed line: approximation H = Go; black solid line: local approximation taking into account 
the boundary. 

Figs. [15] and [16] show the efficiency of this approximation in two different cases: in a 
2D discrete domain, and in a 3D continuous domain. In both cases the approximation 
improves the basic approximation H = Gq, but when the source is near another boundary, 
a systematic deviation appears, due to the influence of the other boundary. 

The curvature of the boundary may be taken into account by approximating the pseudo- 
Green function by the pseudo-Green function inside a circle ( IA6I) or a sphere ( 1A7I) . or outside 
a circle or a sphere (it can be found in 231] ). 



APPENDIX B: COMPUTATION OF THE HIGHER-ORDER MOMENTS 



Discrete case 



In this part we will compute the higher-order moments of the FPT. To do this, we start 
from an extension of Kac's formula (see Appendix [D]), which is the relation between the 
Laplace transforms of the FRT to the subset E, averaged on E, and of the FPT to this 
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Source position x 




FIG. 16: (color online) 3D Brownian motion: the domain is an eight of sphere; the sphere is of 
radius 25, centered on (0,0,0), and the domain is reduced to positive coordinates. The target is 
in (10, 10,2), and the source is in (x,x,3). Red crosses: numerical simulations; blue dashed curve: 
basic approximation H = Gq; solid black curve: approximation (|C7p with H taking into account 
the nearest boundary (Eq. ()A13j) ) 

subset, the starting point being averaged over the complementary subset S. 

tt(E) «e- T > 2 - c -) = (1 - tt(E)) ( C - - l) (e^ (Bl) 
Both averages are weighted by the stationary probability n, in the following sense: 

1 oo 
^ ' iGE i=l 
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-. oo 

< V(T) > 2 = !_ , E x I>fo)X>( T = ( B3 ) 

^ ' i^S t=i 

where Pi(T = t) is the probability for the FRT (or the FPT, according to whether the point 
i belongs to E or not) to be t, if the random walk starts from the point i. To apply the 
equation ( IB II) to the determination of the FPTs, we may notice that the FPT from any 
point of the graph (except target) is the same on the original graph and on the modified 
graph: indeed, the behaviour of a random walk is exactly the same on both lattices as long 
as they do not reach T, and what happens afterwards does not matter. Moreover, the FRT 
to T is still the FPT from S to T, plus one. Thus if we apply the formula (IBlj) to the 
modified graph, E being reduced to T, we get the following relation between the Laplace 
transform of the FPT from S and the FPT averaged over the whole set of points (without 
T): 

J (( e ~ ST >s - 1) = (1 - J) (1 - e s ) <e- sT ) g (B4) 

J is still 7t(iy). We have to pay attention to one thing: the average over E is weighted by 
the weights for the stationary distribution of the modified lattice. To go further we will have 
to consider all the modified lattices with T as target point, the starting point being any 
point of the set. We will denote 7T; the stationary distribution associated with the modified 
graph whose starting point is i, and Jj = 7r;(r T ). Thus, we may note: 

J i « e " ST >„ -!) = (!- eS ) E *'( r i) ( e ~ ST >; ( B5 ) 
From this, we may deduce the recurrence equation for the moments: 



1 m=l ,VT ^ ' 



(B6) 



We may thus compute explicitly the second moment. 



3+T 

If we replace ir and < T > by their values, we get 



' Ji + JiHiTAn) - JiHivAvT)) (H(r T \r T ) - i?(r T | ri )) - ' ' 7 ' 



AT ■ v"V-,,-,y "V-,,-,/, 

(B8) 
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We then may use the value of ^y^, which we know: 

(T 2 ) t = 2N (H(r T \r T ) - H(r T \ri) + Hfaln) - fo|r T )) (H(r T \r T ) - H(r T \ rj )) 

-N{H{t t \t t )-H{t t \t$) (B9) 

This equation is exact, but it is difficult to evaluate properly in the general case. We will 
thus use approximations to evaluate this expression in the case of a 3-D regular lattice, with 
N large and the boundaries far from the target, at a typical distance N^. We can thus 
neglect the term ]V(/7(iy|iy) — i7(r^|rj)) in the R.H.S. of equation (1B9|) . If we develop the 
rest of the formula, we get: 

(T 2 ) t = 2N (NH 2 (r T \r T ) - OT(r T |r r )#(r T | ri ) + F(r T |r T ) ^( H ( r j\ r i) - 2H{r T \r j )) 
+H(r i \r T ) HMtj) - £ if^lr^folr,) + ^ H 2 (r T \r 3 ) ) (BIO) 

We can now drop the least important terms in this formula by evaluating the order of 
magnitude of the various sums over j. We have (cf. Eq. (J6])): 



-Y^H(T i \r j ) = H (Bll) 

j 

Since Gq(t) ~ 1/r in 3D, and the corrections are, on the worst case,of the same order of 
magnitude, we can see that H scales as iV~ 1//3 . If we consider the sums Y2j H 2 ( r T\ r j) an d 
Y2j H( r T\ r j)H(rj\ri), we may first notice that: 

^(rrlriWrjIr,) < #Vl*i) £ i? 2 (r,|r i ) > ) (B12) 

We thus only need to consider the case of (1/N) H 2 (Ti\Tj). And, for the same reasons 
as above, we can see that it scales as N~ 2 ^ 3 . Putting all this together, we have: 

(T 2 ) t = 2N 2 [(H(r T \r T ) - #(r T |r 4 )) (H(r T \r T ) - g) + 0(N- 2 ' 3 )] (B13) 

It is possible to generalize this expression to higher-order moments; we will obtain the 
following result, for a given n: 

(r>, = n\N n [(#(r r |r r ) - H(r T \ ri )) (H(r T \r T ) - H)^ 1 + 0(iV- 2/3 )] (B14) 
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We can prove this by recurrence: if this is true for m < n, then: 



'! I>te) <t»->> 

<Ji 



(B15) 



3+T 



The others terms are negligible (their relative order of magnitude is at most 1 / N) , and we 
will thus ignore them. We replace everything by its value, which gives: 



( 



H(y t \y t ) - R{v T \v>) \ I H(r T \r T 
+H(r j \r i ) - tffolrr) j I -#(r r | rj 




(H(r T \r T ) - H) 
+0(N-l) 



n-2 



(B16) 



Using exactly the same approximations as above (the computation is identical), we get: 



(T n ). = n\N n (H(Y T \v T )-H(v T \v i )){H(Y T \v T )-H) n 1 + 0(N 



-2/3>i 



(B17) 



As for the dependence with n of the correction, since we perform exactly the same operation 
at each step n — > n + 1, the correction will be proportional to n, which may help estimate 
the validity of the approximation. 

This computation fails for elongated domains: two main hypotheses are not satisfied in 
this case, namely that the boundaries are at a typical distance iV 1 / 3 , and that the corrections 
to Go have the same order of magnitude. The method can not either be applied to the 2D 
case, since the terms if 2 (rj|rj) are no longer negligible. 

2. Continuous case 



In the continuous case we can perform a similar computation. The higher-order moments 
of the FPT at the target satisfy the following equations 28]: 



DA(T"(r)) = — n(T (r)) if r G V* 

(T n (r)> = if r G S abs 
<9„(T"(r)> = 0if rGS rcfl 



(B18) 

(B19) 
(B20) 



Using a new time the Green function defined by Eqs. (I36|37II38I) and the Green formula, we 
have 



n 



(T"(r 5 )) = - ^G(r|r s )(T"-' (r)»^r 
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With the knowledge of G(r|r') for all starting points r, it is possible to compute the full 
distribution. In three dimensions, it is possible to find an expression for (T n ) similar to the 
one found in the discrete case. We will start from Eq. fl49l) . We can now compute the second 
moment, using the values for (T) and po : 

(T 2 (r 5 )> = — J ^ [G? (o) + H*(r T \r T ) - H(r T \r s ) + H(r\r s ) - H(r\r T )] (B22) 
[G (a) + H*(r T \r T ) - H(r\r T )} d d v 

To compute this, we will use the two following equations, equivalent to Eq. fIBlip and 
(IB12p for discrete random walks: 

/ H(r \r)d d v = VH + (a 2 ) (B23) 
Jv* 

[ iZ"(r 1 |r)iZ"(r 2 |r)d <i r = O (V 1/3 ) (B24) 
Jv* 

This gives: 

OT/2 

(T 2 (r 5 )) = — [(G (a) + H*(r T \r T ) - H(r T \r s )) (G (a) + H*(r T \r T ) - R) + O (y- 2/3 )] 

(B25) 

This result may be extended by recurrence to higher-order moments, in exactly the same 

way that in the discrete case, which gives: 

nW n r 1 
( T "( r s)> = ~~jy>~ (Go(a) + H *( rT ^ ~ H ( r T\r s )) (G (a) + H*(r T \r T ) - H) n + O («y- 2/3 a 2 " n ) 



(B26) 



APPENDIX C: REFINEMENTS OF THE CONTINUOUS THEORY 

In this Appendix we will see how to improve the results of Section IIII[ provided we know 
the pseudo-Green function H . The results we obtained in Section II III are not perfectly 
satisfying for three reasons: 

• When the source and the target are close, the approximation works better than one 
could naively expect, given that it does not satisfy Eq. (|37l) very accurately. It would 
be interesting to understand why. 

• The approximation lacks accuracy when the target is near a boundary. 
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• In the two-target case the accuracy is not very good when the two targets are close. 

We will treat the first point in detail, and give the corrections, and the method used to 
compute them, for the second and third point. 



1. A better evaluation of G 

To understand this, we will first notice that the Green function we use could also be used 
in an electrostatic problem: the source is equivalent to a point charge, and the absorbing 
spheres are equivalent to conducting spheres set at a null potential. We can thus apply the 



well-known method of images 30| to our problem. If we have an image charge q 



in 3D 



q(r s )={ R " (CI) 

-1 in 2D 

placed on i(rg), located on the line between the center of the sphere and the source, at a 
distance R' = a 2 /R of the target, where R is the source-target distance, then the solution: 

G(r|r s ) = Po (r s ) + G (r|r 5 ) - G (r|r T ) + g(r 5 )(G (r|i(r s )) - G„(r|r T )) (C2) 

satisfies exactly the boundary condition ( |37l) on the target sphere: we have, for r e E a b s 

G (r|r 5 ) - G (r|r r ) + g(r 5 )(G (r|i(r 5 )) - G (r|r T )) = G (r 5 |r T ) - G (a) (C3) 

However, this solution does not satisfy the reflecting boundary conditions, and we will rather 
use the solution: 



G(r|r 5 ) = Po {r s ) + H{r\r s ) - H(t\t t ) + g(r s )(#(r|i(r s )) - H(r\r T )) 



(C4) 



which approximately satisfies fl37|) . provided that we neglect the variations of H*(r\rg) and 
H*(t\tt) on the target sphere. With this approximation we get: 

Po(r 5 ) = G (a) - H(r\r s ) + H*(r T \r T ) + q(r s )(H*(i(r s )\r T ) - H*(r T \r T )) (C5) 

Note that the last term q(rs)(H* (i(r s)\tt)— H*(rT\rT)) can be neglected, since the variations 
of H* over the target sphere are neglected. Finally, to find the Eq. ( )49|) . the only condition is 
to neglect the variations of H* over the target sphere, which will be a good approximation as 
soon as the target is far from any boundary. If this condition is satisfied, the approximation 
for the MFPT is accurate, even if the source is near the target. 
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FIG. 17: (color online) Picture of the real and image charges when the target is near the boundary(+ 
= red plusses; - = blue crosses) 

2. Influence of a boundary 

If the target is near a boundary, however, H* can no longer be considered as constant 
over the target sphere. To have a good approximation of H, one has to decompose the 
function one step further: 

H(t\t') = G (r|r') + G Q (r|s(r')) + H**(r\r'), (C6) 

where s(r) is the point symmetrical to r with respect to the boundary. This simply explicits 
the image charges due to the boundary, which themselves have images on the target sphere. 
The real and image charges are depicted in Fig. [171 

If we take into account all these charges, it is possible to obtain the following expression 
for the MFPT, valid as long as the target sphere does not touch the boundary: : 

(T(r s )> = ^ (G? (o) - H(r T \r s ) + H*(r T \r T ) - K(r s ) - K(s(r s )) + K(s(r T ))) , (C7) 
where K(r) = q(r)(H*{i(r)\r T ) - H*(r T \r T )). 
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3. Two close targets 



The two-target case can be treated likewise: by considering the images of T\ and T 2 on 
the other sphere, it is possible to compute corrections to the terms H i, H 02 , Hi s and Hi 2 
used in Eqs. and (1541) . These correction are: 

H u = H(r Tl \r s ) + g 2 (r 5 )(F(r Tl |i 2 (r 5 )) - H(r Tl \r T2 )), (C8) 

Hqi = Go(a x ) + iy*(r T Jr T J + g 2 (r T J(iy(r Tl |i 2 (r Tl )) - iy(r Tl |r T2 )), (C9) 

# 12 = #(r Tl |r T2 ), (CIO) 

and similar corrections for H 02 and H 2 $. <lk( r ) an d ifc( r ) denote the value and the position 
of the image charge of r inside T^. 

APPENDIX D: PROOF OF KAC'S FORMULA AND OF ITS EXTENSION 
1. The model 

We use the notations of Section [Til R is an arbitrary finite set of points 1, 2, . . . , N, with 
positions rj., r 2 , . . . , r^- is the transition probability from j to i, and we assume that any 
couple of points % and j in R can be joined by at least one succession of links with non-zero 
transition probabilities. 

Among the points of R, we now arbitrarily define a subset E, and note the complementary 
subset E. Practically, the following properties will mostly be interesting if the number of 
points in E is much smaller than the total number N of points, but it is not necessary for 
the definitions. 

With the definitions, the Perron-Frobenius theorem |l| assures that there exists a sta- 
tionary probability vr(rj), which satisfies: 

7r( r i) = Yl u, « 7r ( r i) ( D1 ) 

j£R 

From now on, we will consider that E is absorbing, which means that the particle is absorbed 
as soon as it goes to the subset. However, it may start from it and go away on the following 
step without being absorbed. Thus, we state that, on any state rj, the particle has a 
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probability Pdiji) to be absorbed on its next step equal to: 

Pd(rj) = ^wji (D2) 

2. Obtention of the formula 

Now, the probability p{r i) t) that the conditional particle is adsorbed exactly at time t, 
starting from state i at time 0, obeys the backward equation: 

Pi r i,t) = E^'* _ l ) w ji (D3) 

if t > 2, and 

p(v l ,l)=p d (v i ) (D4) 
As a result, the Laplace transform p of p(rj,t) satisfies: 

p(r i9 s) - e~ s ^m,-; = e" s E^i' s)w jh (D5) 

where p(r,, 1) has been replaced by its value. We multiply this equation by the stationary 
probability vr(rj) and sum up over all i G R. We notice that, from (IDlj) 

J~] WjiTT^j) = 7r(rj) (D6) 

We thus obtain: 

E^^M^-e^E 71 ^') = e- s ^p(r j , S )7r(r J ) (D7) 

We now define two kinds of average for a quantity ip(t): 
(i) the volume average 

-. oo 

(v(T)> £ = E E *) (D8) 

(m) the surface average 

1 oo 

(V(T)> E = -T^yE E ^(*)P('i. *) (D9) 

^ ' ies t=i 
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where tt(S) and vr(S) are the respective stationary probabilities of the volume and the 
surface: 

tt(S) = ;>>(!-,) (D10) 
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(E) = £>(!:«), (Dll) 

and T denotes the absorption time, which corresponds to the FPT to S, or the FRT to 
E, depending on whether the starting point is on E or E. We thus simply get from (1D7|) the 
following equation: 



7T 



(E) <e-* T ) s + tt(E) <e- sT ) g - tt(E) = e^S) (e^g (D12) 



or 



which is the extended Kac's formula, relating the Laplace transforms of the FRTs and the 
FPTs. Thus, for the first moment of T we obtain the very simple and general result: 

(T) E = rrU (D14) 



7T 



(Kac's formula 22]) 



APPENDIX E: SIMULATION METHODS 



1. Random walks 



For random walks we use a method based on the exact enumeration method 311] • The 
exact enumeration method allows one to compute the exact distribution probability up to 
a given time: at each time step (t > 0), we compute the full probability distribution of the 
random walker, using the master equation: 

p(r,t) = - Kr',*-1) (El) 

r'eiV(r) 

p here is the probability of the random walker to be at position r at time t and to never 
have reached the target site before. N(r) is the ensemble of neighbours of r, which includes 
r itself if r is a boundary site. The initial condition is of course p(r, 0) = 5(r, rs). Note that 
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if we set T = S the algorithm will compute the distribution of the FRT. After this first step, 
we have the probability distribution p(t) of the FPT: 

p(t)=p(r T ,t) (E2) 

The last step of the algorithm is to set p(ry,i) to 0, and we can then proceed to the 
computation for the time t+1. This enables us to compute the exact probability distribution, 
but of course the algorithm has to stop at a certain time. To go further, we can notice that the 
tail of the probability distribution is exponential (this corresponds to the highest eigenvalue 
of the transition matrix, the transition probabilities to and out of the target being set to 
to take the absorption into account). If p ~ e~ at for high enough t, then we can compute 
the distribution up to a time to, then estimate: 

<T> = 'gpM + f^ + !^ (E3) 

The two latter terms correspond to Y2t^t p(to)e~ a ^~ to \ Since a is small, its order of mag- 
nitude being l/N, they are approximated by p(t )t /a + p(t ) / a 2 . To estimate a, we take 

1 , p(t - 10) 
a = — In . . E4 
10 pfo) V ' 

(we took 10 steps and not one in order to avoid parity effects). To select to, we run a few 

trial simulations, with a large maximum time, and we determine the minimal to which gives 

a result differing by at most 0.1% from the result obtained with a larger to- We add a small 

security margin, and then run the simulation. We use similar methods for all the other 

quantities studied. The error on the simulation results is thus guaranteed to be less than 

0.1%! 

2. Brownian motion 

Unfortunately, for the Brownian motion, we do not have such an accurate algorithm, and 
we thus used a Brownian-dynamics-based algorithm 32J: we average the time needed to 
reach the target on n = 10 5 Brownian processes. To simulate the Brownian motion, we use 
the following algorithm: 

1. Find the distance between the particle and the nearest obstacle (target, non-flat bound- 
ary). 
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2. Multiply that distance by a constant a (we used a = 0.2) to get a trial typical step 
length. 

3. If we are very close to a boundary, or very close to the target, this trial step length 
would be too small. We thus add a lower cutoff to this trial step length (we took 0.01 
near the target, of typical size radius 1, and 0.2 near the curved boundaries, whose 
radius of curvature was typically 25), and get the typical step length r step . 

4. We use this step length to determine our time step t step = r 2 step . (we have D = 1). 

5. For each direction x,y,z, we add to the position a Gaussian random variable, of 
variance 2t s t ep - To get such a variable, we use two random variables v and /i uniformly 
distributed between and 1, and then the random variable r step ^/ — 2 ln(i/) cos(27r [l) is 
indeed a Gaussian with the required variance. 

6. If we are outside the domain, we move the particle inside the domain, to a position 
symmetrical with respect to the boundary. 

7. If we are inside the target, we end the process, otherwise we take another step. 

This algorithm is less accurate than the one we used in the discrete case, and is compu- 
tationally more expensive. Moreover, the study of the probability density of the FPT is 
delicate with this algorithm. 



APPENDIX F: PROPERTIES OF THE PSEUDO-GREEN FUNCTION H 



The properties of the continuous Pseudo-Green function are well described in 23], and 
we will just describe the properties of the discrete one. We consider the case of symmetric 
transition probabilities. We define the discrete Laplacian operator: 

= - Wij (Fl) 

This operator is hermitian, which will be useful. We define $ p and X p the eigenvectors and 
(real) eigenvalues of the operator —A, ordered from to N — 1 in increasing order. We have 
Aq = 0, and $o — 1/ y/~N, with the usual normalization. Since the operator is hermitian, we 
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can take <3>* = $ p . We define: 



N-l 



(F2) 



This solution satisfies: 



AH(Ti\Tj) = 5., 



1 



(F3) 



which corresponds to the definition we used for H, and we thus found the solution (up to 
a constant) to the equation (jHJ) we used to define H. This shows that H is symmetric in 
its arguments if W = {wij} is symmetric. To prove that the sum Hj = J2iLi -^( r «l r j) i s 
independent of j, we will simply sum up the equation (jSJ) over all i, and use the fact that 
H is symmetric. This gives: 



-AHi 



3 ~ 







(F4) 



and H is proportional to $ 0; an d thus is a constant. 
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